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Abstract 



OO I Because optical systems have huge bandwidth and are capable of generating low noise short pulses they are ideal 



for undersampling multi-band signals that are located within a very broad frequency range. In this paper we propose 

a new scheme for reconstructing multi-band signals that occupy a small part of a given broad frequency range under 

the constraint of a small number of sampling channels. The scheme, which we call multi-rate sampling (MRS), 

entails gathering samples at several different rates whose sum is significantly lower than the Nyquist sampling 

rate. The number of channels does not depend on any characteristics of a signal. In order to be implemented with 

OO I simplified hardware, the reconstruction method does not rely on the synchronization between different sampling 

channels. Also, because the method does not solve a system of linear equations, it avoids one source of lack 

C/j I of robustness of previously published undersampling schemes. Our simulations indicate that our MRS scheme is 

O . robust both to different signal types and to relatively high noise levels. The scheme can be implemented easily 

-^ I with optical sampling systems. 



I. Introduction 



t^ ■ A multi-band signal is one whose energy in the frequency domain is contained in the finite union of 

^ . closed intervals. A sparse signal is a signal that occupies only a small portion of a given frequency region. 

CL(i In many applications of radars and communications systems it is desirable to reconstruct a multi-band 

sparse signal from its samples. When the signal bands are centered at frequencies that are high compared 

to their widths, it is not cost effective and often it is not feasible to sample at the Nyquist rate Fnyq; 

^ the rate that for a real signal is equal to twice the maximum frequency of the given region in which the 

fSj ; signal spectrum is located. It is therefore desirable to reconstruct the signal by undersampling; that is 

C^I to say, from samples taken at rates significantly lower than the Nyquist rate. Sampling at any constant 

rate that is lower than the Nyquist rate results in down-conversion of all signal bands to a low frequency 

region called a baseband. This creates two problems in the reconstruction of the signal. The first is a loss 

56 '. of knowledge of the actual signal frequencies. The second is the possibility of aliasing; i.e. spectrum at 

O ■ different frequencies being down-converted to the same frequency in the baseband. 

>• ■ Optical systems are capable of very high performance undersampling ^. They can handle signals 

K> , whose carrier frequency can be very high, on the order of 40 GHz, and signals with a dynamic range as 

5h ' high as 70 dB. The size, the weight, and the power consumption of optical systems make them ideal for 

- - ■ undersampling. The simultaneous sampling of a signal at different time offsets or at different rates can 

be performed efficiently by using techniques based on wavelength-division multiplexing (WDM) that are 

used in optical communication systems. 

There is a vast literature on reconstructing signals from undersampled data. Landau proved that, 
regardless of the sampling scheme, it is impossible to reconstruct a signal of spectral measure A with 
samples taken at an average rate less than A ^. This rate A is commonly referred to as the Landau rate. 
Much work has been done to develop schemes that can reconstruct signals at sampling rates close to the 
Landau rate. Most are a form of a periodic nonuniform sampling (PNS) scheme [I3]|-|l9l. Such a scheme 
was introduced by Kohlenberg [3] who applied it to a single-band signal whose carrier frequency is known 
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a priori. The PNS scheme was later extended to reconstruct multi-band signals with carrier frequencies 
that are known a priori ||4l, |f8l. 

In a PNS scheme m low-rate cosets are chosen out of L cosets of samples obtained from time- uniformly 
distributed samples taken at a rate F where F is greater or equal to the Nyquist rate Fnyq flU. Consequently, 
the sampling rate of each sampling channel is L times lower than F and the overall sampling rate is L/m 
times lower than F. The samples obtained from the sampling channels are offset by an integral multiple 
of a constant time increment, 1/F. This sampling scheme may resolve aliasing. In a PNS scheme the 
signal is reconstructed by solving a system of linear equations [i4J. PNS schemes can often achieve perfect 
reconstructions from samples taken at a rate that approaches the Landau rate under the assumption that 
the carrier frequencies are known a priori. However, in order to attain a perfect reconstruction, the number 
of sampling channels must be sufficiently high such that the equations have a unique solution flU. 

When the carrier frequencies of the signals are not known a priori, in a PNS scheme, a perfect 
reconstruction requires the sampling rate to exceed twice the Landau rate [|5l, O. In addition, in a 
PNS scheme the number of sampling channels must be sufficiently high BH. Under these two conditions, 
a solution to the set of equations in PNS scheme may be obtained assuming that the sampled signal is 
sparse H. When a PNS scheme is applied to an A^-band real signal (A^ bands in the interval [0, F„y^/2]), 
at least AN channels are required for a perfect reconstruction O, H. A method for obtaining a perfect 
reconstruction has been demonstrated only with the number of channels equal to 8A^ O. Even when the 
requirement of perfect reconstruction is relaxed, the number of channels required to obtain an acceptably 
small error in the reconstructed signal may be prohibitively large. Furthermore, the implementation of 
the schemes to attain the minimum sampling rate relies heavily on assumed values of the widths of the 
sample bands and the number of bands of the signal BH. In the case that the bands of the signal have 
different widths, a PNS scheme for obtaining the minimum sampling rate has not been demonstrated. 

Other important drawbacks of PNS schemes stem from the fact that the systems of equations to be solved 
are poorly conditioned [|7l- Thus, the schemes are sensitive to the bit number of A/D conversion. They 
are also sensitive to any noise present in a signal and to the spectrum of the signal at any frequencies 
outside of strictly defined bands. Moreover, the use of undersampling significantly increases the noise 
in each sampling channel since the noise in the entire sampled spectrum is downconverted to low 
frequencies. Therefore, the dynamic range of of the overall system is limited. The noise may be reduced by 
increasing the sampling rate in each channel. However, since the number of channels needed for a perfect 
reconstruction is determined only by the number of signal bands, the overall sampling rate dramatically 
increases. Another important drawback of PNS scheme is the requirement of a very low time jitter between 
the samplings in the different channels. 

In this paper we propose a different scheme for reconstructing sparse multi-band signals. The scheme, 
which we call multi-rate sampling (MRS), entails gathering samples at P different rates. The number P 
is small (three in our simulations) and does not depend on any characteristics of a signal. Our approach is 
not intended to obtain the minimum sampling rate. Rather, it is intended to reconstruct signals accurately 
with a very high probability at an overall sampling rate that is significantly lower than the Nyquist rate 
under the constraint of a small number of channels. 

The success of our MRS scheme relies on the assumption that sampled signals are sparse. For a typical 
sparse signal, most of the sampled spectrum is unaliased in at least one of the P channels. This is in 
contrast to the situation that prevails with PNS schemes. In PNS schemes, because all channels are sampled 
at the same frequency, an alias in one channel is equivalent to an alias in all channels. 

In our MRS scheme, the sampling rate of each channel is chosen to be approximately equal to the 
maximum sampling rate allowed by cost and technology. Consequently, in most applications, the sampling 
rate is significantly higher than twice the maximum width of the signal bands as usually assumed in PNS 
schemes. 

Sampling at higher rates has a fundamental advantage if signals are contaminated by noise. The spectrum 
evaluated at a baseband frequency /^ in a channel sampling at a rate F is the sum of the spectrum of the 
original signal at all frequencies /^ + mF that are located in the system bandwidth, where m ranges over 



all integers. Thus, the larger the value of F, the fewer terms contribute to this sum. As a result, sampling 
at a higher rate increases the signal to noise ratio in the base-band region. 

To simplify the hardware needed for the sampling, our reconstruction method was developed to not 
require synchronization between different sampling channels. Therefore, our method enables a significant 
reduction in the complexity of the hardware. Moreover, unsynchronized sampling relaxes the stringent 
requirement in PNS schemes of a very small timing jitter in the sampling time of the channels. We also 
do not need to solve a linear set of equations. This eliminates one source of lack of robustness of PNS 
schemes. Our simulations indicate that MRS schemes are robust both to different signal types and to 
relatively high noise. The ability of our MRS scheme to reconstruct parts of the signal spectrum that 
alias when sampled at all P sampling rates can be enhanced by using more complicated hardware that 
synchronizes all of the sampling channels. 

The paper is organized as follows. In section 2 we present some general mathematical background. In 
section 3 we describe the algorithm. In section 4 we give some considerations regarding our algorithm 
complexity. In section 5 we present results of computer simulations. 

II. Mathematical Background and Notation 
A multi-band signal is one whose energy in the frequency domain is contained in a finite union of 
closed intervals IJn=i['^«' ^']- ^ multi-band signal x{t) is said to be sparse in the interval [-Fmin, -flnax] if 
the Lebesgue measure of its spectral support A(x) = J2n=ii^ri — a„) satisfies A <^ F^ax — -^min- 

The signals we consider are sparse multi-band with spectral measure A. We use the following form of 
the Fourier transform of a signal x(t): 

/oo 
xit)expi-2mft). (1) 

-oo 

If the signal x{t) is real (as is every physical signal), then its spectrum X satisfies X{f) = X{—f) where 
a + hi = a — hi and a and h are real numbers. Thus, a real multi-band signal x{t) has fourier transform 
X{f) which, when decomposed into its support intervals, can be represented by 

N 

X{f)=Y,[Sn{f)+^n{-f)\, (2) 

n=l 

where S'„(/) ^ only for / e [a„, 6„] (where 6„ > a„ > 0), and [a„, 6„] n[«m, &m] = for all n ^ m. 

We assume that Fnyq is known a priori. That is to say, we assume that each 6„ for a real signal is at 
most some known value Fnyq/2. Sampling a signal x{t) at a uniform rate F* produces a sampled signal 

oo 

x\t)=x{t + ^^) E '^(^-|l)' (3) 

n=— oo 

where A' is a time offset between the clock of the sampling system and a hypothetical clock that defines 
an absolute time for the signal. Because we are assuming a lack of synchronization between more than one 
sampling channel, we assume that the time offsets A* are unknown. Reconstructing the amplitude of the 
signal spectrum with our scheme does not require knowledge of the time offsets. Only in reconstructing 
the phase of the signal in the frequency domain, do we need in some cases to extract the differences 
between time offsets. 

The Fourier transform of a sampled signal x*(t), X^f), is given by 

oo 

X\f) =F' Y^ ^(f + ^^0 exp[27rz(/ + nF')A% (4) 

n=— oo 

The connection between the spectrum of a sparse signal X(f) and the spectrum of its sampled signal 
X\f) is illustrated in Fig. [B 
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Fig. 1. Illustration of the spectrum of a sparse one-band real signal (a), and the spectrum of its samples which are obtained for the sampling 
rates F^ (b) and F (c). At /o, the signal is unaliased at the sampling rate F^, but is aliased at the sampling rate F . 



One immediate consequence of Eq. |4] is that, up to a phase factor that does not depend on the signal, 
exp[27ri(/ + nF')A*], X*(/) is periodic of period F\ It is also clear that, for a real signal x{t), X*(— /) = 
X^{f). Thus, all of the information about |X*(/)| is contained in the interval [0,F'^/2]. Beside a linear 
chirp caused by the time offset A* all the information about the phase of X\f) is also contained in 
the interval [0,FY2]. We shall refer to this interval [0, FY2] as the i\h baseband. The down-conversion 
of a frequency / G [0, -Fnyq/2] to this baseband is represented by the down-conversion function Z^* : 
[0,F„,q/2]^[0,FV2]: 



D\f) = min[/ mod F\ (F' - /) mod F' 



(5) 



In the case of band-limited signal X(f), for a given frequency /, all but a finite number of terms in the 
infinite sum on the right side of Eq. |4] vanish. If the number of non- vanishing terms is greater than one 
for a given sampling rate F\ then the signal is said to be aliased at / when sampled at the rate F\ If at 
a frequency / only a single term in the sum is not equal to zero, the signal X{f) is said to be unaliased 
at a sampling rate F\ Illustration of aliasing can be seen in Fig. [IJc). In the case of sparse signals, x(t) 
is unaliased at considerable part of it spectral support. The success of an MRS scheme lies in the fact 
that whereas a signal may be aliased at a frequency / when sampled at a rate F\ the same signal may 
be unaliased at the same frequency / when sampled at a different rate F^ . 

Each support interval [a, b] (b > a > 0) of the multi-band signal will be referred to as an originating 
band. According to Eq. H sampling at the rate F' down-converts each originating band [a, b] to a single 
band in the baseband [«*,/?*]. We shall refer to the interval [a*,/?*] as a down-converted band. 

It is apparent that when a single down-converted band [a*,/5*] is given, it is in general not possible to 
identify its corresponding originating band. However, it follows easily from Eq. |4] that the corresponding 
originating band must reside within the set of bands defined by 



Q [a' + mF\ (3' + mF'] ) IJ ( U h^' + ^^'' ""^ + ^^T ) f fl t^' ^"y^/^l ' (6) 




where m is an integer. The set in Eq. [6]can be represented as a finite number of disjointed closed intervals, 
which we denote by [a^, 6^ . We shall refer to each of these intervals as an up-converted band. For clarity, 
we denote all down-converted intervals with greek letters superscripted by the sampling frequency and 
denote all up-converted intervals with latin letters. 

In general, the number of possible originating bands is reduced by sampling at more than one rate. For 
each sampling rate rate F\ an originating band [a, h] must reside within the union of the upconverted 
bands: [a, 6] e U„[a^,6^]. Since the union of upconverted bands is different for each sampling rate, 
sampling at several different rates gives more restrictions over the originating band [a, h]. When sampling 
at P rates, F\ . . . , F^, the originating band must reside within nflj^ U„ [a^, 6^]. 

III. Reconstruction Method 

In this section we describe an algorithm to reconstruct signals from an MRS scheme. First, we describe 
an algorithm for reconstructing ideal multi-band signals, as defined above. Then we present modifications 
to enable a reconstruction of signals that may be contaminated by noise outside of strictly defined bands. 
While such signals are not exactly multi-band, we still consider them multi-band signals provided that 
the noise amplitude is considerably lower than the signal amplitude. 

The reconstruction is performed sequentially. In the first step sets of intervals in the band [0, Fnyq/2] that 
could be the support of X{f) are identified. These are sets that, when down-converted at each sampling 
rate F*, give energy in intervals in the baseband where significant energy is observed. For each hypothetical 
support, the algorithm determines the subsets of the support that are unaliased in each channel. According 
to Eq. m for the correct support, the amplitude of each sampled signal spectrum is proportional to the 
original signal spectrum over the unaliased subset of the support. As a result, for each pair of channels, 
the amplitudes of the two sampled signal spectra are proportional to one another over the subsets of the 
hypothetical support which are unaliased in both channels. Thus, we define an objective function that 
quantifies the consistency between the different channels over mutually unaliased subsets of the support. 
The algorithm chooses the hypothetical support that maximizes the objective function. The amplitude is 
reconstructed from the sampled data on the unaliased subsets of the chosen hypothetical support. In the 
last step, the phase of the spectrum of the originating signal is determined from the unaliased subset of 
the chosen hypothetical support. 

A. Noiseless signals 

In this subsection we assume that all signals are ideal multi-band signals. Although what follows applies 
to more general signals, we assume that all signals have piece-wise continuous spectrum. 

1) Reconstruction of the spectrum amplitude: For each sampled signal X^{f), we consider the indicator 
function X*(/) that indicates over which frequency intervals the energy of the sampled signal X\f) resides. 
To ignore isolated points discontinuity we define the indicator functions X*(/) as follows: 

y , , f 1 for all / e [0, F„yq/2] such that for all e > 0, jj^^ \X%f)\^df' > 
1^ otherwise. 

For piece-wise continuous function, it is simple to show that X'(/) = 1 on closed intervals. 
We define the function X(/) as follows: 

p 
lif) = l[rif), fe[0,F,yj2]. (7) 

Thus, the function X(/) equals 1 over the intersection of all the up-converted bands of the P sampled 
signals. We denote the intervals over which X(/) = 1 by f/i . . . Uk- The Appendix gives sufficient 
conditions under which each originating band coincides with one of the intervals Ui...Uk- Thus, it 
remains to determine which of the K intervals coincide with the originating intervals. 
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For each A; = l,2,---,ifwe consider the indicator function 

1 if / G f7 
otherwise. 
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It follows immediately from Eq. [8] that 



K 



X{f) = J2Uf)- (9) 



fc=l 



To find which sets of Uk (or Ik{f)) match the originating bands each indicator function Xfc(/) is down- 
converted to the baseband via the formula 

(n=oo \ 

Y, Uf + nF^)+U-f + nF')]. (10) 

n=—OQ J 

In Eq. [IO]X[o,FV2](/) is the indicator function of the closed interval [0, FY2]: 

^[o,FV2]UJ-<^ otherwise. ^^^^ 



if (/) is the Heaviside step function 



m) H ? in > 0. <'^> 



The Heaviside step function in Eq. [TO]is used to assure that X^(/) is an indicator function. In the case in 
which the down-conversions of an interval f/^ are aliased at some frequency / within the baseband the 
argument of the step function is an integer greater than 1. However, X^(/) = 1. If, for a frequency / in 
the baseband there is no signal in any of its replicas; i.e., F(nF* ± /) = for all n, then ii(/) = 0. As 
a consequence, X^(/) = also. Therefore, the function X^(/) is equal to one over the down-conversion 
of the interval U^ corresponding sampling rate F*. 

We consider the power set of f/, V{U\, i.e., the set of all subsets of {f/i, ■ ■ ■ , Vk). We denote an 
element of V^} hy U = {f/fci, ■ ■ ■ , Ukg] (0 < Q < K). A subset U G V{U} is deemed to be a support 
consistent combination if, for each sampling rate F\ the down conversion of its intervals matches the 
down-converted bands of the corresponding sampled signal. In terms of indicator functions, we define for 
each U G V{U} the indicator functions 

ruU)=Y.^l^f) /e[o,i^V2]. (13) 

Uk&U 

The function X^(/) is an indicator function for the down-conversion of the intervals of U. Next, we 
define the objective function 

Em = f2 r '\rM)-r{f)\df. (14) 

i=l -^0 

Support consistent combinations are those U for which EiiU) = 0. 

Figure [2] illustrates our method for the signal shown in Fig. [U The support of the signal at positive 
frequencies, shown in Fig. [U consists of a single interval. Figures |2a) and I2b) are graphs of X^(/) and 
X^(/). FigurelUc) is a graph of X(/). The function X(/) is equal to one over four intervals Ui, . . . , U4. Each 
combination of these four intervals must be checked for support consistency. In the example illustrated in 
Fig. I2I we check whether the subset U = {U2} G V{U} is support consistent. Figures |2] (d) and (e) show 
the indicator functions for the down-conversion of U2 at rates F^ and F^: Tij^{f) and Xl,^{f), respectively. 
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Fig. 2. Illustration demonstrating how support-consistency is checked. The input of the algorithm is the sampled signals whose spectra 
X (/) and X (/) are shown Figs. 1 (b) and (c), respectively; their respective indicator functions T^ (/) and T (/) are shown in Fig. 2 (a) and 
(b). Figure 2 (c) shows the indicator function X(/) — X^{f)X^{f). In Figs. 2 (d) and (e), we check whether the subset U = {U2} £ ^{U^ 
is support consistent. Figures [2] (d)and (e) show the indicator functions for the down-conversion of U2 at rates F^ and F^: 1-}j (f) and 
ifj^{f), respectively. The dashed lines illustrate U2, —U2 and their down-conversions. It is evident that the functions X^(/) and Xu^{f) are 
not equal. Hence, U = {U2} i& not a support-consistent combination. 



The dashed lines illustrate U2, —U2, and their down-conversions. It is evident that the functions X^(/) 
and ll/^{f) are not equal. Hence, U = {U2} is not a support-consistent combination. 

Amongst all support consistent combinations U, it is necessary to identify the one that exactly matches 
the originating bands. For this purpose, we introduce two additional objective functions. The support 
consistent combination U that optimizes these function is deemed to be the correct one. 

Amongst support-consistent combinations, amplitude consistent combinations are defined by the am- 
plitudes of the sampled signals at unaliased intervals. Let U = {f/ji, ■ ■ ■ , Uj^} be a support consistent 
combination. Denote the union of all intervals in |jr=i ^jn ^^^^ ^^^ unaliased when down-converted at 
rate F* by S^^ C IJ™=i Uj„. For the correct choice of W, at a frequency / that is unaliased when sampled 
at rates F'^ and F'^ ( f e S^^ n S^^), the functions \X'^{f)\/F'^ and \X'^{f)\/F'^ must be equal. 
Accordingly, we define a second objective function: 



E2{U) =J2[ {\X'^{f)\/F'^ - \X''{f)\/F^n' df. 



(15) 



ji^ia "^w' '^u 



For the correct U, the objective function E2{U) must equal zero. A support-consistent combination U for 
which E2{U) = is said to be amplitude consistent. 



Unfortunately, there may be more than one amplitude-consistent combination. This is the case, for 
example, when for all zi and 12, S^ fl S^ is empty. In such cases, the objective function E2{U) cannot 
be sufficient to identify the correct U. Thus, we introduce a third objective function E^{U). This function 
favors options in which the integrals in Eq. [15] are calculated over large sets. The third objective function 
is defined by 

E3(W)=5^A(S-nS-), (16) 

where A(S[J fl S^) is the Lebesgue measure of S^ fl S^. The amplitude-consistent combination that 
maximizes E-i[U) is deemed to be the correct one. In the rare case that E^iU) is maximized by more 
than one amplitude-consistent combination, the outcome of the algorithm is not determined. 

After the optimal U = {f/ji, ■ " " ) ^jm} is chosen, the amplitude of the signal is reconstructed from 
the samples. We define the function r(/) as the number of sampled signals which are unaliased at the 
frequency /: r(/) = J2i=i'^T.^ (/)' where I-^t (/) is the indicator function of the interval S^, defined 
similarly to Eq. \TT\ For each / within the detected originating bands, i.e. / G ljr=i ^in' if '^(/) > 0' we 
reconstruct the corresponding amplitude of the spectrum at / from the sampled signals by 

1 ^\X%f)\I^.{f) 

In words, for each frequency / that is unaliased in at least one channel, the signal amplitude is averaged 
over all the channels that are not aliased at /. For all other frequencies, notably those that alias in all 
sampling channels, Xu{f) is set to equal zero. 

2) Reconstruction of the spectrum phase: The spectrum of a signal can be expressed as X{f) = 
|X(/)| exp{j arg[X(/)]}. In the previous section we described how to reconstruct the amplitude |X(/)| 
from the signal's sampled data. In this section we describe a method of reconstructing the phase arg[X(/)]. 
If the time offsets A* of Eq. |4] were known a priori, reconstructing the phase would be trivial. The 
reconstruction in this case could be performed by using a variant of Eq. \T7} with |X*"(/)| replaced 
by X*"(/) exp(— 27r/A*"). This would yield a full reconstruction of the signal (phase and amplitude). 
However, because of the lack of synchronization between the channels, the time offsets A* are not known 
a priori. Consequently, it is more difficult to reconstruct the phase. After identifying the signal bands, we 
can calculate the differences A*i — A*^ between two different time offsets. This is sufficient to enable the 
reconstruction of the phase of the signal spectrum up to a single linear phase factor. 

The difference between two time offsets A*^ and A*^ can be calculated directly in the case that S[J fl S^ 
contains at least one finite interval. In this interval the phase of X*^(/)/X*2(/) satisfies the following 
equation: 

arg[X*^(/)/X^^(/)] = 2nf{A'' - A'^) + 27rk, for some integer k (18) 

The left side of Eq. [T8] is determined by the sampled data. By performing a linear fit we calculate the 
difference between the two offsets A*^ and A*^. We do this for all pairs of offsets for which S^ fl S^ 
contains at least one finite interval. 

There may exist cases in which there exist zi and Z2 such that S^ fl S^ does not contain one finite 
interval but for which A*^ — A*^ can still be calculated. For example, in the case of three offsets A*^ , A*^ 
and A*^', if one can calculate (A'^ — A*^) and (A*^ — A*^), then (A*^ — A*^) can also be calculated by 
simple algebra. If there exist in, ■ ■ ■ im, such that for each n < k < m. — 1, Yl'-^ fl S^^^ contains at least 
one finite interval, then we say that in and im are phase connected and denote this by z„ ~ z^- If ^ ~ j? 
then difference between the two offsets A-^ — A* can be calculated. In the case S^^ does not contain any 
finite intervals, we define A* ~ A*. It is clear that ~ is an equivalence relation ifTOl and thus partitions 
the A* into equivalence classes. 

For each A'^ and A*^ in the same class, one can calculate their difference. One can obtain a full 
reconstruction of the phase if there exists one class C such that each originating frequency is unaliased in 



at least one channel belonging to C; i.e, there exist a class C = A*" . . . A*"\ such that UfcLn ^u = Ufc=i ^jk^ 
where U = {Uj^, ■ ■ ■ , Uj^}. 

B. Physical signals 

To sample realistic signals (i.e., not strictly multi-band and in the presence of noise), the algorithm 
needs to be adjusted. In this subsection we describe adjustments to our algorithm to overcome the noise. 
The algorithm requires five new parameters. In section 5, we give examples of reconstructing signals 
contaminated by strong noise. In those examples, the success of the reconstruction does not depend on 
the exact choice of the five parameters. 

In the presence of noise, the definition of the support of the sampled signals must be adjusted. First, 
a small ^ is chosen. Then, a small positive threshold value T is chosen. The indicator function X*(/) is 
then redefined as follows: 

j.(j) ^ I 1 if / G [0, F„,q/2] and ^ //_+/ {X^'M' > T ^^^^ 

\ otherwise. 

The choice of the threshold T depends on the average noise level. 

When reconstructing physical signals, it is not reasonable to expect EiiU) to equal for any combination 
U. An initial adjustment is to require that Ei{U) < b for some positive b. The shortcoming of this condition 
is that the threshold b does not depend on the signal. To make the threshold to depend on the signal in a 
simple way, we introduce the following condition: 

Ei{U) < a mm [Ei{U)] + b (20) 

where a > 1 is a chosen parameter. The parameters a and b control the tradeoff between the chance 
of success and runtime. If a and b are too small, the correct subset U may not be included in the set 
of support constituent combinations. On the other hand, if a and b are too large, then the number of 
support-consistent combinations may be large. This results in a slow run time. 

Finally, we make two modifications to the objective function E^{U). We replace the length of the 
mutually unaliased intervals by a weighted energy of the sampled signals in these interval. The objective 
function E^(V() is replaced with E^: 

2 x^^(f) ' 



Em = E / 

• /• -'0 



pii 



W,,^i,{f,U)df, (21) 



where IVj^ J2(/, W) is a weight function. The weight function favors combinations in which the sampled 
signals are similar in mutually unaliased internals and is defined in the following. 

We first note that for each two channels ii and i2, the intersection of their non-aliased supports (S^ flS^) 
is a union of a finite number of disjoint intervals V^/^'^^ ■ ■ ■ V^'^^. We define 

_ /^n.,||X-^(/)|/F-^-|X-^(/)|/F'Md/ 

^"'"^^^ " ' /^H,.||Xn(/)| + |XM/)||d/ • ^^^^ 

Finally, we define the weight function: 

W^n..(/) = J2^M-Pl^Lnm^v:^'^^if)^ (23) 

k 

where p is a chosen positive constant and Tyii.i2{f) is the indicator function of the interval V^^'''^. The 
parameter p is chosen according to an assumed signal to noise ratio (SNR). When the SNR is lower, in 
order to accept higher errors p is chosen to be smaller. In the case of a noiseless signal and an amplitude- 
consistent U, each /i*^^ ^^ vanishes. Therefore, in this case, each element in the sum on the right-hand side 
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of Eq. [2T| gives the energy of the signal over S^ fl S^. In all other cases, the energy in each interval 
V;!^'*^ is weighted according to the relative error between X*^(/) and X^^{f) over V^!^'*^ 

Since in the case of noisy signals, neither Ei(U) nor E2{U) vanishes for the combination which 
corresponds to the originating bands, both Ei(U) and E2(U) should be considered in the final step of 
choosing the best combinations. Accordingly, we define the following objective function E^oiiU): 

mmu {El iU) } mmu {^2 iU) } ^^^^ I ^^ (^) 



for all U such that mmu {Ei{U)},mmu{E2iU)}-,"Ci\iiiu [E^iiU)} 7^ 0. Amongst all such combinations 
that also satisfy Eq. |20l the one that gives the maximum value of EioilU) is deemed to be correct. In cases 
in which either minj/ {EiiU)}, mmu {E2{U)} or min^^ {E^{U)} equals zero for a certain combination hi, 
the maximum of E-i(U) is chosen as the solution. 

To reconstruct the phase, the only change made is in how the difference between the offsets is calculated. 
Equation [TSl holds for all the disjoint intervals V^*^'*^ G S]^ fl S^. Accordingly, we perform the linear 
fit for each intervals, and obtain a certain value for A'^ — A'^ . Each value is weighted by the length of 
its respective V^!^'*^. These weighted values are averaged. The result is an estimate for A'^ — A*^. This 
averaging procedure may increase the accuracy in the estimate of A*i — A*^ . 

IV. Complexity considerations 

In this section we discuss considerations used to reduce the computational complexity of our algorithm. 
Choosing a subset U E ViU} involves calculating three objective functions. We explain why eliminating 
possibilities through the use of EiilA) alone can significantly reduce runtime. 

In the first step of the algorithm, we find support consistent combinations by calculating the objective 
function Ei{U) for elements in V{U}. Assuming the largest element in V{U} contains K intervals, and 
that the signal is composed of up to N bands in [0,Fnyq/2], the number of elements in V{U} that one 
needs to check is equal to 

n=l ^ ^ 



In the case N ^ K, the complexity is approximately 0(2^). When N/K <^ 1, the last term in Eq. 
number of options to be checked is approximately equal to 0{K^ /N\). 

The complexity of checking a single option out of V{U} for support consistency (Eq. [HI) is 0(1) and 
it does not depend on the number of points used to discretize the spectrum. By contrast, the complexity 
of checking such an option for amplitude consistency (Eqs. [15] and \T6^ is of the order of the number of 
points used to represent the spectrum. This is a major reason for using the support-consistency criterion 
to narrow down the number of options needed to be checked for amplitude consistency. The amplitude 
consistency is calculated only for support-consistent options, which are in general much fewer than what 
is prescribed by Eq. 



V. Numerical Results 

This section describes results of our numerical simulations. The simulations were carried out in the 
two cases considered in the previous sections: i) ideal multi-band signals and ii) noisy signals. In all our 
examples, the number of channels P was set equal to three, P = 3. 

In all our simulations, the number the bands in [0, Fayq/2] equals N, where A^ < 4. Unless stated 
otherwise the band number refers to the number of bands in the non negative frequency region [0, Fnyq/2]. 
Using the notations in Eq. [21 each signal in each band is given by 

^ , , f A„ cos[7r(/ - fn)/Bn] if 2|/ - /„|/B„ < 1 
"■^•'^ 10 otherwise, 
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where Bn is the spectral width of the nth band, /„ is its central frequency, and An is the maximum 
amplitude. The total spectral measure of the signal support equals S^. = 2 J2n=i ^"' ^^^ ^^^ minimal 
sampling rate is equal to 2Sa, [[6], twice the Landau rate. In each simulation, all the bands had the same 
width, i.e. i?„ = T.x/{2N). The amplitudes An were chosen independently from a uniform distribution 
on [1,1.2]. The central frequencies /„ were also chosen independently from a uniform distribution on the 
region [0,Fnyq/2]. We eliminated cases in which there was an overlap between any two different bands. 
The time offsets. A* were chosen independently from a uniform distribution on [0, l/B„]. 

In each of the simulations, we set B = 800 MHz and 40 < F^yq < 76 GHz. This choice of parameters is 
consistent with previous optical sampling experiments [IJ. The sampling rates were chosen as F^ = 3.8Fo, 
i^2 _ 4^^^ ^^^ p3 _ 4.2Fo, where the value Fq varied between simulations. These sampling rates were 
chosen such that, for each pair of sampling rates (F*, F^), the functions X*(/), T^{f) do not have a 
common multiple smaller than Fnyq. This condition is satisfied for all Fq > Fnyq/76. 

To obtain an exact reconstruction, the resolution in which the spectrum is represented A/ should be 
such that the discretization of the originating baseband downconverts exactly to the discretization grid in 
each baseband. This condition is satisfied when F^Af (i = 1,2,3) is an integer. In our examples, we 
used a spectral resolution A/ = 0.8MHz for all the channels. 

The use of the same spectral resolution for all channels is not only convenient for implementation of 
our algorithm, but it also compatible with the operation of the sampling system used in our experiments 
[[ll. In the implementation of the sampling system, an optical system performs the down conversion of 
the signal by multiplying it by a train of short optical pulses. In each channel a different repetition rate 
of the optical pulse train is used. The sampled signal in each channel is then converted into an electronic 
signal and passed through a low-pass filter which rejects all frequencies outside the baseband. The P 
filtered sampled signals have a limited bandwidth. These signals are sampled once more, this time at a 
constant rate, using P electronic analog to digital converters. The use of the optical system allows the use 
of electronic analog to digital converters whose bandwidth is significantly lower than the bandwidth of 
the multi-band signal [IJ. Because the signals at the basebands are sampled with the same time resolution 
and have the same number of samples, their spectra, which are obtained using the Fast Fourier Transform, 
have the same spectral resolution. 

In the first set of simulations we increased the signal bandwidth, without changing the sampling rates. 
We used two performance criteria: correct detection of the originating bands and exact reconstruction of the 
signal. As to the first criterion, we required only that the spectral support of the signal be detected without 
an error. As to the second criterion, we required that the signal spectrum (phase and amplitude) be fully 
and exactly reconstructed without any error. Because the second criterion concerns exact reconstructions, 
in the case that the algorithm failed to reconstruct the signal at even a single frequency, it was considered 
to have failed the second criterion. 

We chose Fq = 1 GHz. This corresponds to a total sampling rate Ftot = F^ + F'^ + F^ which equals 
15 times the Landau rate (7.5 the minimum possible rate). The statistics were obtained by averaging over 
1000 runs. Figures [3] (a) and (b) show the results for signals with 3 and 4 positive bands, respectively, 
as a function of the Nyquist rate. In Fig. [3] (a), the percentage of a correct band detection is shown by 
the squares, whereas the full reconstruction percentage is shown by circles. The open circles and squares 
represent the results obtained when the maximum number of bands assumed by the algorithm was 3, 
and the dark circles and squares represent the cases in which the maximum assumed band number was 
equal to 4. Figure |3](b) shows the band-detection percentage (solid curve) and reconstruction percentages 
(dashed curve) in the case that both the maximum number of originating and assumed bands is 4. The 
figures show that both the success percentages were high and were not significantly dependent on the 
Nyquist rate of the signal or on the number of assumed bands. 

Figure |4] shows the average run time as a function of the Nyquist rate. The results in the case of 4 input 
bands in which the assumed maximum number of bands is 4 is shown in the solid curve. The results in 
the case of 3 input bands is shown with the dotted curve in the case of 3 assumed bands and with the 
dashed curve in the case of 4 assumed bands. The results show that while an increase in the Nyquist rate 
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Fig. 3. Success percentage for the first set of simulations with Fq = 1 GHz as a function of Nyquist rate. In Fig. [3] (a), the percentage 
of a correct band detection is shown by the squares. The full reconstruction percentage is shown by circles. The open circles and squares 
represent the results obtained when the assumed maximum number of positive bands equals 3. The dark circles and squares represent the 
cases in which the maximum assumed positive band number equals 4. Figure [5] (b) shows the band-detection percentage (solid curve) and 
reconstruction percentages (dashed curve) in the case that both the maximum number of originating and assumed positive bands equals 4. 
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Fig. 4. The run time for the second set of simulations as a function of the Nyquist rate. The results in the case of 4 input positive bands 
with assumed number of positive bands equals 4 is shown in the solid curve. The results in the case of 3 input positive bands is shown with 
the dotted curve in the case of 3 assumed positive bands and with the dashed curve in the case of 4 assumed positive bands. 



does not significantly affect the reconstruction statistics, it results in an increase in run time. 

In the second set of simulations, we measured the performance of our algorithm as a function of Fq. 
The Nyquist rate used in the simulation was F^yq = 40 GHz. For each choice of Fq, the statistics were 
obtained by averaging over 500 runs. The results did not change significantly when the averaging was 
performed over 1000 runs. The simulation was run for the same number of originating bands and assumed 
bands as in the first set of simulations. Figures [5] (a) and (b) show the success percentages for signals 
with 3 and 4 bands, respectively, and Fig. [6] shows the average run time. The two success percentages and 
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run time are shown as a function of the total sampling rate Ftot divided by the Landau rate, FLandau = 800 
MHz. The symbols used in Figs. [5] (a) and (b) and Fig. [6] correspond to those used in Figs. |3](a) and (b) 
and Fig. Irrespectively. 
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Fig. 5. Success percentage for the first set of simulations as a function of the sum of the sampling rates divided by the Landau rate. As 
in Fig. [3] in Fig. [5] (a), the percentage of a correct band detection is shown by the squares. The full reconstruction percentage is shown by 
circles. The open circles and squares represent the results obtained when the assumed maximum number of positive bands equals 3. The dark 
circles and squares represent the cases in which the maximum assumed positive band number equals 4. Figure |5](b) shows the band-detection 
percentage (solid curve) and reconstruction percentages (dashed curve) in the case that both the maximum number of originating and assumed 
positive bands equals 4. \ 



The results shown in Figs. [5] (a) and (b) demonstrate that, in all the cases that we checked, the average 
percentage of successful band detection was over 99.5% for sampling frequencies above 8 times the 
Landau rate. The reconstruction percentages were lower than these band-detection percentages and were 
also much more affected by the sampling rate and by the number of originating bands. As expected, the 
run time increases dramatically with reduction of the sampling rate and also increases with the assumed 
maximum number of bands. We ran similar simulations with different numbers of originating bands and 
different numbers of assumed bands. The trends were similar. 

In the final set of simulations, the signals are noisy. We added to the originating signal white Gaussian 
noise in the band [— Fnyq/2, Fnyq/2], where Fnyq = 40 GHz. We denote by a the standard deviation of the 
Gaussian noise in the pre-sampled signal. Upon sampling the signal at rate F*, the standard deviation of 



the noise increases to a* = crA/fFnyq/F'] owing to aliasing of the noise, where \x\ equals the smallest 
integer greater or equal to x. 

In the this set of simulations, we reconstructed signals with different noise levels added. We chose 
^ = 6 MHz. The threshold was chosen to be T = 2maxj(cr*). Accordingly, the parameter p in Eq. [23] was 
chosen to be p = maxj(a*). The other parameters used in the simulation were a = 2 and 6 = 16 MHz. 
Because the signals were not ideal, an exact reconstruction was not possible and the definitions of an 
accurate band detection and accurate reconstruction needed to be changed. A band detection was deemed 
accurate if the originating bands approximately matched the reconstructed bands. A signal reconstruction 
was deemed accurate if the signal's originating bands were detected accurately and if each reconstructed 
band Xu(f) satisfied 

f |X^(/)-X(/)|<max(aOBm. (27) 

Here X(f) is the noiseless signal and the integration is performed over only the detected band. In a correct 
reconstruction, it is expected that the average reconstruction error is lower than the standard deviation of 
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Fig. 6. The run time for the first set of simulations as a function of the sum of the sampling rates divided by the Landau rate in the cases 
of signals with 4 and 3 positive bands. The results in the case of 4 input positive bands with assumed number of positive bands equals 4 is 
shown in the solid curve. The results in the case of 3 input positive bands is shown with the dotted curve in the case of 3 assumed positive 
bands and with the dashed curve in the case of 4 assumed positive bands. 
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of the added noise. The figure shows the band-detection percentage (solid curve) and reconstruction percentages (dashed curve) in the case 
that both the maximum number of originating and assumed positive bands equals 4. 



the noise in the noisiest channel, i.e. the channel at the lowest sampling rate. We chose the same sampling 
rates as those chosen in the second set of simulations. For these rates: maxj(cr*) = 3.3a. 

The detection percentages and reconstruction percentages are shown in Fig. Ul The figure clearly shows 
that high percentages are obtained even in the case of low signal to noise ratio. We repeated this last 
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set of simulations using Gaussian signals instead of the signals of Eq. |26l We found that results are not 
sensitive to the specific choice of signal type. 

VI. Conclusion 

Typical undersampling schemes are PNS schemes. In such schemes samples are taken from several 
channels at the same low rate. These schemes have many drawbacks. In this paper we propose a new 
scheme for reconstructing multi-band signals under the constraint of a small number of sampling channels. 
We have developed an MRS scheme; a scheme in which each channel samples at a different rate. We 
have demonstrated that sampling with our MRS scheme can overcome many of the difficulties inherent 
in PNS schemes and can effectively reconstruct signals from undersampled data. For a typical sparse 
multi-band signal, our MRS scheme has the advantage over PNS schemes because in almost all cases, 
the signal spectrum is unaliased in at least one of the channels. This is in contrast to PNS schemes. With 
PNS schemes an alias in one channel is equivalent to an alias in all channels. 

Our MRS scheme uses a smaller number of sampling channels than do PNS schemes. We also choose 
to sample at a higher sampling rate than PNS schemes use in attaining the theoretical minimum overall 
sampling rate required for a perfect reconstruction. The use of higher rates has an inherent advantage in 
that it increases the sampled signal to noise ratio. Our MRS scheme also does not require the solving 
of poorly conditioned linear equations that PNS schemes must solve. This eliminates one source of lack 
of robustness of PNS schemes. Our simulations indicate that our MRS scheme, using a small number of 
sampling channels (3 in our simulations) is robust both to different signal types and to relatively noisy 
signals. 

Our reconstruction scheme does not require the synchronization of different sampling channels. This 
significantly reduces the complexity of the sampling hardware. Moreover, asynchronous sampling does 
not require very low jitter between the sampling time at different channels as is required in PNS schemes. 
Our reconstruction scheme resolves aliasing in almost all cases but not all. In rare cases, reconstruction of 
the originating signal fails owing to aliasing. One of the methods to resolve aliasing is to synchronize the 
sampling in all the channels. With such synchronization, aliasing can be resolved by inverting a matrix 
similarly to as is done in PNS schemes. However, such an approach requires both much more complex 
hardware and a larger number of sampling channels that sample with a very low jitter. Moreover, in case 
of signals that are aliased simultaneously in all channels, the noise in the reconstructed signal is expected 
to be much stronger than the noise in the original signal. 

Future work should focus on testing our algorithm's ability to reconstruct experimental data. Optical 
systems for performing experiments are currently in existence. 

VII. Appendix 

In section 3.A.1 we have denoted the intervals over which the indicator function X(/) = Ihy Ui . . . Uk- 
In this appendix we give a sufficient and necessary conditions under which the spectral support of a signal 
coincides with a subset U of {Ui, . . . Uk} and under which the function Ei(U) (Eq. [141) is equal to zero. 
Although it applies for more general cases, we assume that the function X{f) is piecewise continuous. 

The conditions are as follows: 

1) For each frequency /q which fulfills //"J*"^^ |X(/)p(i/ > for all £ > 0, we obtain that 

jJlTe \X\f)\'^df > for all £ > and 1 < z < P. 

2) For each originating band with support [a, h], there exists an interval [a — e, a + e], (e 7^ 0) whose 
down-converted band does not overlap any other down-converted band in at least one of the sampled 
signals. Similarly, for each originating band with support [a, h], there exists an interval \b — e,h + e\, 
whose down-converted band does not overlap any other down-converted band in at least one of the 
sampled signals. 



16 

Condition 1 assures that originating bands are contained within U^^f/j. Condition 2 guarantees that the 
originating coincide exactly with a subset of V{U}. It is obvious that when the conditions are satisfied, 
Ei{U) = 0. 

The first condition excludes cases in which the down-converted bands cancel each other's energy over 
a certain interval due to destructive interference. When the condition is fulfilled, for each frequency /o 
within the originating bands, we obtain X(/o) = 1. Thus, each originating band [a,b] is contained within 
one of the intervals that make up the support of 1(f). Mathematically, for each [a, b], there exist Uk, such 
that [a,b] C Uk- 

The second conditions assures us that for each originating band [a, b], the intervals [a — e, a] and [b, b + e] 
are not contained within any of the Uk for all values of e. Consequentially, if [a,b] C Uk, then [a,b] = Uk- 
When the two conditions are fulfilled, we obtain that there exist a set of intervals U, which matches the 
originating bands, and for which EiiU) = 
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